participants <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/sub_ses_passQC_finallist.csv", header = FALSE)
cifti.vertex.mask <- ciftiTools::read_xifti("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/study_cortexmask_vertex.dscalar.nii")
glasser.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/glasser360_regionlist.csv", header = T)
schaefer.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/schaefer400_regionlist.csv", header = T)
SAaxis <- read.csv("/Users/valeriesydnor/Documents/ResearchProjects/Neuron Neurodevelopmental Imaging/CorticalMap_Code/Sensorimotor_Association_Axis_AverageRanks.csv", header=T)
brains <- read.csv("/Users/valeriesydnor/Documents/ResearchProjects/Neuron Neurodevelopmental Imaging/CorticalMap_Code/brainmaps_glasser.csv", header=T)
mappys <- read.csv("/Users/valeriesydnor/Software/mappys/mappings_400.csv", header=T)

Parcellate REHO Data

Parcellate vertex level REHO data into glasser360 and schaefer400 parcels for analysis

alpraz_reho_parcellate <- function(subid, sesid){
  
  #read in xcp_abcd left and right hemisphere reho giftis
  reho.lh <- readgii(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/xcpabcd/task_regress/surface/xcp_abcd/%1$s/%2$s/func/%1$s_%2$s_task-emotionid_space-fsLR_den-32k_hemi-L_desc-reho_den-91k_bold.func.gii", subid, sesid)) #32492 vertex data in $reho.lh$data$normal
  reho.rh <- readgii(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/xcpabcd/task_regress/surface/xcp_abcd/%1$s/%2$s/func/%1$s_%2$s_task-emotionid_space-fsLR_den-32k_hemi-R_desc-reho_den-91k_bold.func.gii", subid, sesid)) #32492 vertex data in $reho.rh$data$normal
  
  #combine left and right hemisphere data and save as cifti
  xii.reho <- as_xifti(cortexL = reho.lh$data$normal, cortexL_mwall = cifti.vertex.mask$meta$cortex$medial_wall_mask$left, cortexR = reho.rh$data$normal, cortexR_mwall = cifti.vertex.mask$meta$cortex$medial_wall_mask$right)
  dir.create(path = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/", subid, sesid), recursive = TRUE)
  write_cifti(xii.reho, file.path(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_fsLRvertex.dscalar.nii", subid, sesid)))
  
  #parcellate vertex reho data with workbench
  #GLASSER360 HCP Multimodal template parcels
  command1 = sprintf("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_fsLRvertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_glasser360.pscalar.nii", subid, sesid)
  ciftiTools::run_wb_cmd(command1, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
  
  #SCHAEFER400 template parcels
  command2 = sprintf("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_fsLRvertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/Schaefer2018_400Parcels_17Networks_order.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_schaefer400.pscalar.nii", subid, sesid)
  ciftiTools::run_wb_cmd(command2, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
}
for(row in c(1:nrow(participants))){
  subid=participants[row,1]
  sesid=participants[row,2]
  alpraz_reho_parcellate(subid, sesid)
}

Study-Specific Cortical Mask

Create a parcellated study-specific group mask based on the fslr vertex-level group mask

#parcellate the vertex-level group mask 

#GLASSER360
ciftiTools::run_wb_cmd("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/study_cortexmask_vertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_glasser360_unthresholded.pscalar.nii", intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)

#SCHAEFER400
ciftiTools::run_wb_cmd("-cifti-parcellate /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/study_cortexmask_vertex.dscalar.nii /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/Schaefer2018_400Parcels_17Networks_order.dlabel.nii COLUMN /cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_schaefer400_unthresholded.pscalar.nii", intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
#threshold the group mask and save out binary 0/1 mask vector

#GLASSER360
cifti.glasser.mask.unthresholded <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_glasser360_unthresholded.pscalar.nii")
#information about the unthresholded mask: 191 parcels have 100% of vertices included in the group mask, 3 parcels have 50-60% of vertices included, 6 parcels have 60-70%, 5 parcels have 70-80% of vertices included, 9 parcels have 80-90% of vertices included, 18 parcelts have 90-99% of vertices included

#threshold the parcellated mask to retain only parcels where >= 50% of vertices have data for all subjects (final parcel number = 232)
glasser360.parcelmask.unthresholded <- cifti.glasser.mask.unthresholded$data
glasser360.parcelmask.thresholded <- ifelse(glasser360.parcelmask.unthresholded < 0.5, 0, glasser360.parcelmask.unthresholded)
glasser360.parcelmask.thresholded <- ifelse(glasser360.parcelmask.thresholded >= 0.5, 1, glasser360.parcelmask.thresholded)

#SCHAEFER400
cifti.schaefer.mask.unthresholded <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_schaefer400_unthresholded.pscalar.nii")
#information about the unthresholded mask: 178 parcels have 100% of vertices included, 5 parcels have 50-60% of vertices, 2 parcels have 60-70% of vertices, 12 parcels have 70-80%, 9 havae 80=90%

#threshold the parcellated mask to retain only parcels where >= 50% of vertices have data for all subjects (final parcel number = 227)
schaefer400.parcelmask.unthresholded <- cifti.schaefer.mask.unthresholded$data
schaefer400.parcelmask.thresholded <- ifelse(schaefer400.parcelmask.unthresholded < 0.5, 0, schaefer400.parcelmask.unthresholded)
schaefer400.parcelmask.thresholded <- ifelse(schaefer400.parcelmask.thresholded >= 0.5, 1, schaefer400.parcelmask.thresholded)
#save mask as vector
write.csv(x = glasser360.parcelmask.thresholded, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_glasser360_thresholded.txt", row.names = F, quote = F)

write.csv(x = schaefer400.parcelmask.thresholded, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/atlases/study_cortexmask_schaefer400_thresholded.txt", row.names = F, quote = F)

Visualize mask

#visualize the parcellated and thresholded group mask

#GLASSER360
finalmask <- as.data.frame(glasser360.parcelmask.thresholded)
finalmask$label <- glasser.parcel.labels$label
ggseg(.data = finalmask, atlas = "glasser", mapping=aes(fill=glasser360.parcelmask.thresholded), position = c("stacked")) + theme_void() + scale_fill_gradient2(low= "black", high = "dodgerblue3", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.5) + theme(legend.position = "none")

rm(finalmask)

#SCHAEFER400
finalmask <- as.data.frame(schaefer400.parcelmask.thresholded)
finalmask$label <- schaefer.parcel.labels$label
ggseg(.data = finalmask, atlas = "schaefer17", mapping=aes(fill=schaefer400.parcelmask.thresholded), position = c("stacked")) + theme_void() + scale_fill_gradient2(low= "black", high = "mediumpurple3", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.5) + theme(legend.position = "none")

rm(finalmask)

Create Masked and Parcellated REHO Data Matrices

alpraz_reho_glasser <- function(subid, sesid){
  
  #read in parcellated reho data
  cifti.reho.glasser <- read_cifti(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_glasser360.pscalar.nii", subid, sesid))
  reho.glasser.unmasked <- as.array(cifti.reho.glasser$data)
  reho.glasser.masked <- sweep(reho.glasser.unmasked, MARGIN = 1, glasser360.parcelmask.thresholded, `*`) 
  reho.glasser.masked <- as.data.frame(reho.glasser.masked)
  reho.glasser.masked[reho.glasser.masked == 0] <- NA
  write.table(reho.glasser.masked, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_glasser360_masked.csv", subid, sesid), row.names = F, col.names = F, quote = F)
  return(reho.glasser.masked)
}
reho.subxparcel.matrix.glasser <- matrix(data = NA, nrow = 84, ncol = 364)
regionheaders <- glasser.parcel.labels$orig_parcelname
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))
colnames(reho.subxparcel.matrix.glasser) <- colheaders

for(row in c(1:nrow(participants))){
  subid=participants[row,1]
  sesid=participants[row,2]
  drug=participants[row,3]
  meanFD=participants[row,4]
  data <- alpraz_reho_glasser(subid, sesid)
  reho.subxparcel.matrix.glasser[row,] <- cbind(subid, sesid, drug, meanFD, t(data))
}
write.csv(reho.subxparcel.matrix.glasser, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_glasser.csv", row.names = F, quote = F)
alpraz_reho_schaefer <- function(subid, sesid){
  
  #read in parcellated reho data
  cifti.reho.schaefer <- read_cifti(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_schaefer400.pscalar.nii", subid, sesid))
  reho.schaefer.unmasked <- as.array(cifti.reho.schaefer$data)
  reho.schaefer.masked <- sweep(reho.schaefer.unmasked, MARGIN = 1, schaefer400.parcelmask.thresholded, `*`) 
  reho.schaefer.masked <- as.data.frame(reho.schaefer.masked)
  reho.schaefer.masked[reho.schaefer.masked == 0] <- NA
  write.table(reho.schaefer.masked, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/%1$s/%2$s/%1$s_%2$s_reho_schaefer400_masked.csv", subid, sesid), row.names = F, col.names = F, quote = F)
  return(reho.schaefer.masked)
}
reho.subxparcel.matrix.schaefer <- matrix(data = NA, nrow = 84, ncol = 404)
regionheaders <- schaefer.parcel.labels$label
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))
colnames(reho.subxparcel.matrix.schaefer) <- colheaders

for(row in c(1:nrow(participants))){
  subid=participants[row,1]
  sesid=participants[row,2]
  drug=participants[row,3]
  meanFD=participants[row,4]
  data <- alpraz_reho_schaefer(subid, sesid)
  reho.subxparcel.matrix.schaefer[row,] <- cbind(subid, sesid, drug, meanFD, t(data))
}
write.csv(reho.subxparcel.matrix.schaefer, file = "/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_schaefer.csv", row.names = F, quote = F)

Visualize Cortical REHO Maps

Glasser: Average map

Visualze across-cortex regional homogeneity (Glasser)

REHO.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_glasser.csv", header=T)
REHO.glasser.map <- REHO.glasser %>% dplyr::select(-subid, -sesid, -drug, -meanFD)
REHO.glasser.map <- colMeans(REHO.glasser.map)
REHO.glasser.map <- as.data.frame(REHO.glasser.map)
REHO.glasser.map$label <- glasser.parcel.labels$label

ggseg(.data = REHO.glasser.map, atlas = "glasser", mapping=aes(fill=REHO.glasser.map), position = c("stacked")) + theme_void() + scale_fill_gradientn(colors = viridis_pal(option="B")(10), limits = c(.45,.765))

Glasser: S-A axis correlations

Correlation with S-A Axis

cor.test(REHO.glasser.map$REHO.glasser.map[181:360], SAaxis$final.rank, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.glasser.map$REHO.glasser.map[181:360] and SAaxis$final.rank
## S = 276670, p-value = 0.6957
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.03654347
cor.test(REHO.glasser.map$REHO.glasser.map[1:180], SAaxis$final.rank, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.glasser.map$REHO.glasser.map[1:180] and SAaxis$final.rank
## S = 278704, p-value = 0.2896
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.09959796

Correlation with G1 fMRI Ratio

cor.test(REHO.glasser.map$REHO.glasser.map[181:360], brains$G1.fMRI, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.glasser.map$REHO.glasser.map[181:360] and brains$G1.fMRI
## S = 313580, p-value = 0.05948
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1748265
cor.test(REHO.glasser.map$REHO.glasser.map[1:180], brains$G1.fMRI, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.glasser.map$REHO.glasser.map[1:180] and brains$G1.fMRI
## S = 313654, p-value = 0.01075
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2374892

Correlation with T1T2 Ratio

cor.test(REHO.glasser.map$REHO.glasser.map[181:360], brains$T1T2ratio, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.glasser.map$REHO.glasser.map[181:360] and brains$T1T2ratio
## S = 206588, p-value = 0.01443
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2260187
cor.test(REHO.glasser.map$REHO.glasser.map[1:180], brains$T1T2ratio, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.glasser.map$REHO.glasser.map[1:180] and brains$T1T2ratio
## S = 188002, p-value = 0.005444
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2582577

Schaefer: Average map

Visualze across-cortex regional homogeneity (Schaefer)

REHO.schaefer <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho_subxparcel_matrix_schaefer.csv", header=T)
REHO.schaefer.map <- REHO.schaefer %>% dplyr::select(-subid, -sesid, -drug, -meanFD)
REHO.schaefer.map <- colMeans(REHO.schaefer.map)
REHO.schaefer.map <- as.data.frame(REHO.schaefer.map)
REHO.schaefer.map$label <- schaefer.parcel.labels$label

ggseg(.data = REHO.schaefer.map, atlas = "schaefer17", mapping=aes(fill=REHO.schaefer.map), position = c("stacked")) + theme_void() + scale_fill_gradientn(colors = viridis_pal(option="B")(10), limits = c(.45,.765))

Schaefer: S-A axis correlations

Correlation with G1fMRI

cor.test(REHO.schaefer.map$REHO.schaefer.map, mappys$SA_1, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.schaefer.map$REHO.schaefer.map and mappys$SA_1
## S = 2358642, p-value = 0.001502
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2098851

Correlation with T1T2 Ratio

cor.test(REHO.schaefer.map$REHO.schaefer.map, mappys$myelin, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  REHO.schaefer.map$REHO.schaefer.map and mappys$myelin
## S = 1470286, p-value = 0.0001918
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2458045

Placebo versus GABA Agonist Paired T-tests

Glasser

Effect of GABA agonist on regional homogeneity per parcel (Glasser)

P-values

drug <- REHO.glasser %>% filter(drug == 0)
drug <- drug %>% select(-subid, -sesid, -drug, -meanFD)
placebo <- REHO.glasser %>% filter(drug == 1)
placebo <- placebo %>% select(-subid, -sesid, -drug, -meanFD)

#run paired t-tests on each parcel (column)
tstats <- col_t_paired(drug, placebo)
num <- nrow(as.data.frame(tstats$pvalue) %>% filter(`tstats$pvalue` < 0.05))
print(sprintf("Number of significant parcels (uncorrected): %s", num))
## [1] "Number of significant parcels (uncorrected): 42"

P-value Histograms

#p-value histograms
ps.fdr <- p.adjust(na.omit(tstats$pvalue))
par(mfrow=c(1,2))
hist(tstats$pvalue, col = "mediumpurple4", xlab = "Paired t-test p-values", main="")
hist(ps.fdr, col = "dodgerblue4", xlab = "Paired t-test FDR-corrected p-values", main="")

T-statistic Histogram

#distribution of t-values
hist(tstats$statistic, col = "plum2", xlab = "Paired t-test t-values", main="")

T-statistic Map

REHO higher in GABA > placebo, the t-statistic is positive (yellow) REHO lower in GABA < placebo, the t-statistic is negative (purple)

tstatistics <- tstats$statistic
tstatistics <- as.data.frame(tstatistics)
tstatistics$label <- glasser.parcel.labels$label

ggseg(.data = tstatistics, atlas = "glasser",  mapping=aes(fill=tstatistics), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.88)

Correlation with S-A Axis

Correlation with S-A Axis (Glasser)

cor.test(tstatistics$tstatistics[181:360], SAaxis$final.rank, method=c("spearman"))
## Warning in cor.test.default(tstatistics$tstatistics[181:360],
## SAaxis$final.rank, : Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  tstatistics$tstatistics[181:360] and SAaxis$final.rank
## S = 373871, p-value = 7.575e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4007081
cor.test(tstatistics$tstatistics[1:180], SAaxis$final.rank, method=c("spearman"))
## Warning in cor.test.default(tstatistics$tstatistics[1:180], SAaxis$final.rank, :
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  tstatistics$tstatistics[1:180] and SAaxis$final.rank
## S = 378506, p-value = 2.119e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4933579

Correlation with T1T2 Ratio (Glasser)

cor.test(tstatistics$tstatistics[181:360], brains$T1T2ratio, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tstatistics$tstatistics[181:360] and brains$T1T2ratio
## S = 199462, p-value = 0.006101
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2527162
cor.test(tstatistics$tstatistics[1:180], brains$T1T2ratio, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tstatistics$tstatistics[1:180] and brains$T1T2ratio
## S = 171112, p-value = 0.0004238
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3248954

Schaefer

Effect of GABA agonist on regional homogeneity per parcel (Schaefer)

P-values

drug <- REHO.schaefer %>% filter(drug == 0)
drug <- drug %>% select(-subid, -sesid, -drug, -meanFD)
placebo <- REHO.schaefer %>% filter(drug == 1)
placebo <- placebo %>% select(-subid, -sesid, -drug, -meanFD)

#run paired t-tests on each parcel (column)
tstats <- col_t_paired(drug, placebo)
## Warning: col_t_paired: 173 of the columns had less than 2 paired observations.
## First occurrence at column 23
num <- nrow(as.data.frame(tstats$pvalue) %>% filter(`tstats$pvalue` < 0.05))
print(sprintf("Number of significant parcels (uncorrected): %s", num))
## [1] "Number of significant parcels (uncorrected): 54"

P-value Histograms

#p-value histograms
ps.fdr <- p.adjust(na.omit(tstats$pvalue))
par(mfrow=c(1,2))
hist(tstats$pvalue, col = "mediumpurple4", xlab = "Paired t-test p-values", main="")
hist(ps.fdr, col = "dodgerblue4", xlab = "Paired t-test FDR-corrected p-values", main="")

T-statistic Histogram

#distribution of t-values
hist(tstats$statistic, col = "plum2", xlab = "Paired t-test t-values", main="")

T-statistic Map

REHO higher in GABA > placebo, the t-statistic is positive (yellow) REHO lower in GABA < placebo, the t-statistic is negative (purple)

tstatistics <- tstats$statistic
tstatistics <- as.data.frame(tstatistics)
tstatistics$label <- schaefer.parcel.labels$label

ggseg(.data = tstatistics, atlas = "schaefer17",  mapping=aes(fill=tstatistics), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0.88)

Correlation with S-A axis

Correlation with G1 fMRI (Schaefer)

cor.test(tstats$statistic, mappys$SA_1, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tstats$statistic and mappys$SA_1
## S = 2921866, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4987956

Correlation with T1T2 Ratio (Schaefer)

cor.test(tstats$statistic, mappys$myelin, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tstats$statistic and mappys$myelin
## S = 1134078, p-value = 7.194e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4182652

Correlation with PVALB (Schaefer)

cor.test(tstats$statistic, mappys$PVALB, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tstats$statistic and mappys$PVALB
## S = 1305640, p-value = 4.154e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3302611

Correlation with SST (Schaefer)

cor.test(tstats$statistic, mappys$SST, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  tstats$statistic and mappys$SST
## S = 2558606, p-value = 1.776e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.3124583

REHO PCA and Logistic Regressions

Glasser

REHO PCA and PC scores logistic regression for drug/placebo prediction (Glasser)

PCA Variance Explained

inputmatrix <- REHO.glasser
inputmatrix <- inputmatrix %>% dplyr::select(-subid, -sesid, -drug, -meanFD) 
inputmatrix <- inputmatrix %>% mutate_if(is.character,as.numeric)
inputmatrix <- as.data.frame(inputmatrix) #convert matrix to df
inputmatrix[inputmatrix == 0] <- NA
inputmatrix.cleaned <- inputmatrix[, unlist(lapply(inputmatrix, function(x) !all(is.na(x))))] #remove all columns (parcels) with null data (all NAs) #84x232
reho.pca.glasser <- prcomp(inputmatrix.cleaned, scale. = TRUE, center = TRUE) 
print(summary(reho.pca.glasser)$importance[2,1:10]*100)
##    PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8    PC9   PC10 
## 25.324  5.610  5.445  4.826  3.371  3.138  2.807  2.740  2.340  2.259
pve <- 100*(reho.pca.glasser$sdev)^2/sum ((reho.pca.glasser$sdev)^2)
plot(pve, pch=16, 
   xlab="Principal Components",
  ylab="Prop. of variance explained")

Logistic Regression

reho.pca.glasser.scores <- as.data.frame(reho.pca.glasser$x[,1:10])
reho.pca.glasser.scores$subid <- participants$V1 
reho.pca.glasser.scores$sesid <- participants$V2
reho.pca.glasser.scores$drug <- participants$V3 
reho.pca.glasser.scores$meanFD <- participants$V4 #match subids to scores

summary(glm(as.factor(reho.pca.glasser.scores$drug) ~ reho.pca.glasser.scores$PC1 + reho.pca.glasser.scores$meanFD, family=binomial(logit)))
## 
## Call:
## glm(formula = as.factor(reho.pca.glasser.scores$drug) ~ reho.pca.glasser.scores$PC1 + 
##     reho.pca.glasser.scores$meanFD, family = binomial(logit))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.42708  -1.15218   0.06032   1.15661   1.40665  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)
## (Intercept)                     0.65963    0.60666   1.087    0.277
## reho.pca.glasser.scores$PC1     0.02048    0.02952   0.694    0.488
## reho.pca.glasser.scores$meanFD -3.30648    2.84188  -1.163    0.245
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.45  on 83  degrees of freedom
## Residual deviance: 114.26  on 81  degrees of freedom
## AIC: 120.26
## 
## Number of Fisher Scoring iterations: 4

Schaefer

REHO PCA and PC scores logistic regression for drug/placebo prediction (Schaefer)

PCA Variance Explained

inputmatrix <- REHO.schaefer
inputmatrix <- inputmatrix %>% dplyr::select(-subid, -sesid, -drug, -meanFD) 
inputmatrix <- inputmatrix %>% mutate_if(is.character,as.numeric)
inputmatrix <- as.data.frame(inputmatrix) #convert matrix to df
inputmatrix[inputmatrix == 0] <- NA
inputmatrix.cleaned <- inputmatrix[, unlist(lapply(inputmatrix, function(x) !all(is.na(x))))] #remove all columns (parcels) with null data (all NAs) #84x232
reho.pca.schaefer <- prcomp(inputmatrix.cleaned, scale. = TRUE, center = TRUE) 
print(summary(reho.pca.schaefer)$importance[2,1:10]*100)
##    PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8    PC9   PC10 
## 26.054  6.463  5.366  4.711  3.922  3.170  2.893  2.654  2.543  2.226
pve <- 100*(reho.pca.schaefer$sdev)^2/sum ((reho.pca.schaefer$sdev)^2)
plot(pve, pch=16, 
   xlab="Principal Components",
  ylab="Prop. of variance explained")

Logistic Regression

reho.pca.schaefer.scores <- as.data.frame(reho.pca.schaefer$x[,1:10])
reho.pca.schaefer.scores$subid <- participants$V1 
reho.pca.schaefer.scores$sesid <- participants$V2
reho.pca.schaefer.scores$drug <- participants$V3 
reho.pca.schaefer.scores$meanFD <- participants$V4 #match subids to scores

summary(glm(as.factor(reho.pca.schaefer.scores$drug) ~ reho.pca.schaefer.scores$PC1 + reho.pca.schaefer.scores$meanFD, family=binomial(logit)))
## 
## Call:
## glm(formula = as.factor(reho.pca.schaefer.scores$drug) ~ reho.pca.schaefer.scores$PC1 + 
##     reho.pca.schaefer.scores$meanFD, family = binomial(logit))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.46086  -1.14744   0.04421   1.14899   1.38829  
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)                      0.63966    0.60842   1.051    0.293
## reho.pca.schaefer.scores$PC1     0.02505    0.02963   0.845    0.398
## reho.pca.schaefer.scores$meanFD -3.20496    2.84927  -1.125    0.261
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.45  on 83  degrees of freedom
## Residual deviance: 114.02  on 81  degrees of freedom
## AIC: 120.02
## 
## Number of Fisher Scoring iterations: 4

SVM

Functions

SVM_2class <- function(df,folds,feature_selection = F,feature_proportion = .1,num_repetitions = 100){
  #SVM 2-class classifier
  # folds: number of folds for cv.
  # feature_selection: optional flag for data-driven feature selection. Selects the proportion of features indicated by feature_proportion. Not used by default.
  # num_repetitions: How many times to you want to repeat the cv process using different random splits of the data. This is just extra caution against a randomly good or bad split.
  
  cat('\nRunning SVM models.....')
  # Set up folds for CV
  if (folds == "LOO") {
    # This is for leave-one-out cross-validation
    num_folds = length(unique(df$subid))
    num_repetitions <- 1
  } else {
    num_folds = folds
    num_repetitions <- num_repetitions
  }

  svm_output <- vector(mode = "list",length = num_repetitions) #set up output object.
  
  # Create the folds
  unique_IDs <- unique(df$subid)
  subid_folds <- replicate(num_repetitions,sample(unique_IDs,size = length(unique_IDs))) # create sets of random draws.
  foldIdxs <- data.frame(subid=unique_IDs)
  foldIdxs$foldID <- row_number(foldIdxs$subid)
  foldIdxs$foldID <- ntile(foldIdxs$foldID,num_folds)
  # cat('Sending data to CV')
  for (r in 1:num_repetitions) {
    foldIdxs$subid <- subid_folds[,r] # Grab a random split for folds.
    fold_output<-vector(mode = "list",length = num_folds)
    cat(sprintf('\nrepetition  %d.... ',r))
    for (fold in 1:num_folds) {
      # cat(sprintf('\nfold  %d.... ',fold))
      trainingIDs <- as.matrix(foldIdxs %>% filter(foldID != fold) %>% select(subid))
      trainingIndex <- df$subid %in% trainingIDs # indices for training subs
      trainingData <- df[trainingIndex, 3:dim(df)[2] ] %>% arrange(drug) #this is important because libsvm automatically makes the first observation class 1, so drug must be first every time. Placebo will be class -1.
      testData <- df[!trainingIndex, 4:dim(df)[2]] # test data. Take columns 4:end (Removes subid sesid drug).
      testLabels <- data.frame(df[!trainingIndex,c(1:3) ])
      # svm
      x <- as.matrix(trainingData[, 2:dim(trainingData)[2]])
      y <- as.factor(as.matrix(trainingData[,1]))
      svm.model <- svm(x =x, y = y, cost = 1, kernel = "linear",type = "C-classification",scale = F)
      svm.pred <- predict(svm.model, as.matrix(testData))
      
      w <- t(svm.model$coefs) %*% svm.model$SV #calculate feature weights.
      # num_features <- dim(x)[2]
      decisionValues <- w %*% t(testData)-svm.model$rho # Get decision valus.
      distance <- decisionValues/norm(w) #calculate distance from the classification hyperplane
      # just adding the results to the dataframe.
      testLabels$decisionValues <- t(decisionValues)
      testLabels$distance <- t(distance)
      testLabels$model.pred = svm.pred
      fold_output[[fold]]<-testLabels
    }
    svm_output[[r]] <- data.table::rbindlist(fold_output) # saving output for the repetition.
    cat('complete\n')
  }
  # output_list <- apply(subid_folds,CV_function,df=df,num_folds = num_folds,MARGIN = 2)
  
  # Now train a model using all the data. This is for the estimation of the feature weights using the most possible data.
  final_data<-df[, 3:dim(df)[2] ]
  x <- as.matrix(final_data[, 2:dim(final_data)[2]])
  y <- as.factor(as.matrix(final_data[,1]))
  svm.model <- svm(x = x, y = y, 
                   cost = 1, kernel = "linear",type = "C-classification",scale = F)
  w <- t(svm.model$coefs) %*% svm.model$SV
  svm_results <- list(svm_output,w,svm.model)
  cat('\nFinished SVM models\n')
  
  return(svm_results)
}
AUC <- function(DecisionValues, labels){
  # This function creates an ROC curve and then calculates AUC.
  # Decision values is nx1
  # Labels is nx1
  
  # N.B.
  # Drug (class 0) is assigned as +1 by libsvm and placebo (class 1) is assigned as 0 by libsvm default. 
  # Adjust the true labels and predicted labels to match that here so that the decision values make sense.
  labels <- labels*-1+1
  
  P <- sum(labels == 1)
  N <- sum(labels == 0)
  Sorted_DecisionValues <- sort(unique(DecisionValues), decreasing = FALSE)
  numDecisionValues <- length(Sorted_DecisionValues)
  
  TP_Array <- vector(mode = "numeric",length = numDecisionValues)
  FP_Array <- vector(mode = "numeric",length = numDecisionValues)
  Accuracy_Array = vector(mode = "numeric",length = numDecisionValues)
  for (i in 1:numDecisionValues){
    thisCutoff <- Sorted_DecisionValues[i]
    thisPredictedLabels <- as.numeric(DecisionValues>thisCutoff)
    detections <- thisPredictedLabels==1
    
    TP <- sum(labels[detections] == thisPredictedLabels[detections])
    TPR <- TP/P
    FP <- sum(labels[detections]!=thisPredictedLabels[detections])
    FPR <- FP/N
    
    TP_Array[i] <- TPR
    FP_Array[i] <- FPR
    
    Accuracy_Array[i] = (TP + N - FP) / (P + N)
  }
  
  ROC_output <- data.frame(TPR =TP_Array,FPR=FP_Array,Accuracy = Accuracy_Array)
  ROC_output <- ROC_output%>%arrange(TPR,FPR)
  
  #AUC
  dFPR <- c(0,diff(ROC_output$FPR))
  dTPR <- c(0,diff(ROC_output$TPR))
  AUC <- sum(ROC_output$TPR * dFPR) + sum(dTPR * dFPR)/2
  return(AUC)
}

Glasser

Drug v. Placebo SVM (Glasser)

regionheaders <- glasser.parcel.labels$label
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))

df <- REHO.glasser
colnames(df) <- colheaders
df <- df %>% dplyr::select(-meanFD)
df <- df[, unlist(lapply(df, function(x) !all(is.na(x))))]
model_results <- SVM_2class(df, "LOO") # model_results contains [[1]] svm_output table of subid, sesid, drug, decisionValues, distance, and model.pred, [[2]] weights for each region from the model run with all sessions and [[3]] the model itself

saveRDS(model_results, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.glasser.rds")
model_results <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.glasser.rds")
prediction_output <- model_results[[1]]
W <- model_results[[2]] #weights for each parcel
svm.model <- model_results[[3]] #the model run 
num_features <- sum(!is.na(W[1,]))

Model accuracy

accuracy_fun <- function(x) sum(x$model.pred==x$drug)/dim(x)[1] # function to calculate the accuracy. 
accuracies <- sapply(prediction_output,accuracy_fun) #Apply the function to each repetition of the cross-validation. For LOO there is only 1 repetition, and accuracies = accuracy
accuracy <- mean(accuracies)
accuracy
## [1] 0.6190476

Number of correct classifications

num_obs <- length(prediction_output[[1]]$model.pred)
num_correct <- round(accuracy*num_obs)
num_correct
## [1] 52

Binomial p-value significance

b <- binom.test(num_correct,num_obs,.5) #This is a binomial test. The p-value is not used (permutation test is used instead), but this function summarizes the data nicely.
b
## 
##  Exact binomial test
## 
## data:  num_correct and num_obs
## number of successes = 52, number of trials = 84, p-value = 0.03753
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5065562 0.7228868
## sample estimates:
## probability of success 
##              0.6190476
b$pred_data <- prediction_output
b$accuracy <- accuracy

Permutation p-value significance

num_permutations = 1000
nw=1
folds="LOO"

perm_acc <- matrix(nrow = num_permutations)
perm_W <- matrix(nrow = num_permutations,ncol = num_features)
perm_list <- list()

perm_list = foreach(perm_chunk = idiv(num_permutations,chunks = nw),
                        .combine=c,
                        .export = c("featureExtraction","SVM_2class"),
                        .packages = c("dplyr","e1071","matrixTests")) %do% {  # This must be `dopar` to be parallel.
                          perm_result_list=list()
                          for (p in 1:perm_chunk){
                            # thisPerm <- perms[p,]
                            cat(sprintf("\nPermutation %d",p))
                            thisDF <- df
                            # thisDF$drug <- df$drug[thisPerm]
                            permuted <- df %>% select(subid,drug) %>% group_by(subid) %>% mutate(perm_drug=sample(drug))
                            thisDF$drug <- permuted$perm_drug
                            perm_pred_result <- SVM_2class(thisDF,folds,
                                                           feature_selection = feature_selection,
                                                           feature_proportion = feature_proportion,
                                                           num_repetitions = 1)
                        
                            perm_result_list[[p]] = list(perm_pred_result)
                          }
                          perm_result_list
                        }

  #organizing the output
pred_data_list <- lapply(perm_list, function(x) x[[1]][[1]])
perm_Ws <- lapply(perm_list,function(x) x[[1]][[2]])
# Calculating permutation accuracies and AUCs
perm_acc_distribution <- sapply(pred_data_list, function(x) sum(x[[1]]$model.pred==x[[1]]$drug)/length(x[[1]]$drug))

perm_p <- sum(perm_acc_distribution>accuracy)/length(perm_acc_distribution)

perm_auc_distribution <- sapply(pred_data_list, function(x) AUC(DecisionValues=x[[1]]$decisionValues,labels=x[[1]]$drug))
    
    
W_test <- sapply(perm_Ws, FUN = function(x) {
    abs(x)>abs(W)
    })
W_sig <- rowMeans(W_test)
    
b$perm_p <- perm_p
b$perm_W_sig <- W_sig
b$perm_accs = perm_acc_distribution
b$perm_aucs = perm_auc_distribution

  
saveRDS(b, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.glasser.rds")
saveRDS(W, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.glasser.rds")
b <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.glasser.rds")
W <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.glasser.rds")
print("How many times is the permuted accuracy greater than the actual accuracy?")
## [1] "How many times is the permuted accuracy greater than the actual accuracy?"
sum(b$perm_accs>accuracy)
## [1] 1
print(sprintf("Overall Accuracy: %1.3f; p = %.5f\n\n",accuracy,b$perm_p))
## [1] "Overall Accuracy: 0.619; p = 0.00100\n\n"

P-value distribution

hist(b$perm_accs, col="grey")
abline(v = accuracy, col="dodgerblue3", lwd=3, lty=2)

Model weights map

SVM.weights.glasser <- as.data.frame(t(W))
SVM.weights.glasser$label <- row.names(SVM.weights.glasser)
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_p9.46v"] <- "rh_R_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_a9.46v"] <- "rh_R_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_9.46d"] <- "rh_R_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_OP2.3"] <- "rh_R_OP2-3"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_p9.46v"] <- "lh_L_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_a9.46v"] <- "lh_L_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_9.46d"] <- "lh_L_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_OP2.3"] <- "lh_L_OP2-3"


ggseg(.data = SVM.weights.glasser, atlas = "glasser",  mapping=aes(fill=V1), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL)

Model weights S-A axis correlation

SAaxis$label <- glasser.parcel.labels$label[181:360]
SVM.weights.glasser.SA <- merge(SVM.weights.glasser, SAaxis, by="label")
cor.test(SVM.weights.glasser.SA$V1, SVM.weights.glasser.SA$final.rank, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.glasser.SA$V1 and SVM.weights.glasser.SA$final.rank
## S = 343225, p-value = 0.00178
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2858925

Model accuracy compared to Catch22 Features

catch22_SVM_accuracy <- function(feature){
  catchnames <- (c("DN_HistogramMode_5","DN_HistogramMode_10","CO_f1ecac","CO_FirstMin_ac",
              "CO_HistogramAMI_even_2_5", "CO_trev_1_num", "MD_hrv_classic_pnn40", 
              "SB_BinaryStats_mean_longstretch1", "SB_TransitionMatrix_3ac_sumdiagcov",
              "PD_PeriodicityWang_th0_01", "CO_Embed2_Dist_tau_d_expfit_meandiff",
              "IN_AutoMutualInfoStats_40_gaussian_fmmi", "FC_LocalSimple_mean1_tauresrat",
              "DN_OutlierInclude_p_001_mdrmd", "DN_OutlierInclude_n_001_mdrmd", 
              "SP_Summaries_welch_rect_area_5_1", "SB_BinaryStats_diff_longstretch0",
              "SB_MotifThree_quantile_hh", "SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1",
              "SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1", "SP_Summaries_welch_rect_centroid",
              "FC_LocalSimple_mean3_stderr"))
  
  #prepare data for SVM
  catchname <- catchnames[feature]
  catchdata <- read.csv(sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/catch%s_sesxparcel_matrix.csv", feature))
  df1 <- catchdata[,1:3]
  df2 <- catchdata[,5:364]
  df2 <- df2 %>% mutate_if(is.character,as.numeric)
  df <- cbind(df1, df2)
  df <- df[, unlist(lapply(df, function(x) !all(is.na(x))))]
  
  #Run SVM
  catchmodel_results <- SVM_2class(df, "LOO") # model_results contains [[1]] svm_output table of subid, sesid, drug, decisionValues, distance, and model.pred, [[2]] weights for each region from the model run with all sessions and [[3]] the model itself
  
  saveRDS(catchmodel_results, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_modelresults_catch%s.Rdata", feature))

  catchprediction_output <- catchmodel_results[[1]]
  catchprediction_weights <- catchmodel_results[[2]] #weights for each parcel

  #Calculate accuracy and binomial p-value
  accuracy_fun <- function(x) sum(x$model.pred==x$drug)/dim(x)[1] # function to calculate the accuracy. 
  accuracy <- sapply(catchprediction_output,accuracy_fun)
  
  num_obs <- length(catchprediction_output[[1]]$model.pred)
  num_correct <- round(accuracy*num_obs)
  b <- binom.test(num_correct,num_obs,.5) #This is a binomial test
  binomial_pvalue <- b$p.value
  
  #Calculate permutation-based p-values
  
  num_permutations = 1000
  nw=1
  folds="LOO"
  num_features = 232

  perm_acc <- matrix(nrow = num_permutations)
  perm_W <- matrix(nrow = num_permutations,ncol = num_features)
  perm_list <- list()

  perm_list = foreach(perm_chunk = idiv(num_permutations,chunks = nw),
                        .combine=c,
                        .export = c("featureExtraction","SVM_2class"),
                        .packages = c("dplyr","e1071","matrixTests")) %do% {  # This must be `dopar` to be parallel.
                          perm_result_list=list()
                          for (p in 1:perm_chunk){
                            # thisPerm <- perms[p,]
                            cat(sprintf("\nPermutation %d",p))
                            thisDF <- df
                            # thisDF$drug <- df$drug[thisPerm]
                            permuted <- df %>% select(subid,drug) %>% group_by(subid) %>% mutate(perm_drug=sample(drug))
                            thisDF$drug <- permuted$perm_drug
                            perm_pred_result <- SVM_2class(thisDF,folds,
                                                           feature_selection = feature_selection,
                                                           feature_proportion = feature_proportion,
                                                           num_repetitions = 1)
                        
                            perm_result_list[[p]] = list(perm_pred_result)
                          }
                          perm_result_list
                        }

  
  pred_data_list <- lapply(perm_list, function(x) x[[1]][[1]])
  perm_Ws <- lapply(perm_list,function(x) x[[1]][[2]])
  
  perm_acc_distribution <- sapply(pred_data_list, function(x) sum(x[[1]]$model.pred==x[[1]]$drug)/length(x[[1]]$drug))
  saveRDS(perm_acc_distribution, file = sprintf("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_permutedaccuracies_catch%s.Rdata", feature))
  
  perm_p <- sum(perm_acc_distribution>accuracy)/length(perm_acc_distribution)
  perm_auc_distribution <- sapply(pred_data_list, function(x) AUC(DecisionValues=x[[1]]$decisionValues,labels=x[[1]]$drug))

  compiled.output <- cbind(catchname, accuracy, binomial_pvalue, perm_p)
  
  return(compiled.output)
}
catch22.SVM.accuracy <- matrix(NA, ncol = 4, nrow = 22)
for(feature in c(1:22)){
  thisfeature.accuracy <- catch22_SVM_accuracy(feature)
  catch22.SVM.accuracy[feature,] <- thisfeature.accuracy
}
colnames(catch22.SVM.accuracy) <- c("Catch Feature", "SVM Accuracy", "Binomial pval", "Permutation pval")
write.csv(catch22.SVM.accuracy, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/Catch22_SVM_Accuracies.csv", quote = F, row.names = F)
catch22.SVM.accuracy
catch22.SVM.accuracy <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/Catch22_SVM_Accuracies.csv")
hist(as.numeric(catch22.SVM.accuracy[,2]), col="grey", breaks=12, xlim = c(0.4,0.65))
abline(v = accuracy, col="red", lwd=3, lty=2)

Best Catch22 Features

Catch9

catch9.modelresults <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_modelresults_catch9.Rdata")
prediction_output <- catch9.modelresults[[1]]
W <- catch9.modelresults[[2]]
SVM.weights.glasser <- as.data.frame(t(W))
SVM.weights.glasser$origlabel <- row.names(SVM.weights.glasser)
SVM.weights.glasser$label <- NA
for(weight in c(1:232)){
  newlabel <- substr(SVM.weights.glasser$origlabel[weight],1,nchar(SVM.weights.glasser$origlabel[weight])-4)
  if(substring(newlabel, 1, 1) == "R"){
    newlabelhemi <- gsub("R_", "rh_R_", newlabel)}
  else
    newlabelhemi <- gsub("L_", "lh_L_", newlabel)
SVM.weights.glasser[weight,3] <- newlabelhemi
}
SVM.weights.glasser.SA <- merge(SVM.weights.glasser, SAaxis, by="label")
cor.test(SVM.weights.glasser.SA$V1, SVM.weights.glasser.SA$final.rank, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.glasser.SA$V1 and SVM.weights.glasser.SA$final.rank
## S = 269943, p-value = 0.1958
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1225927
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_p9.46v"] <- "rh_R_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_a9.46v"] <- "rh_R_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_9.46d"] <- "rh_R_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_OP2.3"] <- "rh_R_OP2-3"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_p9.46v"] <- "lh_L_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_a9.46v"] <- "lh_L_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_9.46d"] <- "lh_L_9.46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_OP2.3"] <- "lh_L_OP2-3"
ggseg(.data = SVM.weights.glasser, atlas = "glasser",  mapping=aes(fill=V1), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL)
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  side  region label           V1 origlabel  
##   <chr> <chr> <chr> <chr> <chr>  <chr>        <dbl> <chr>      
## 1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_9.46d -0.0154 L_9.46d_ROI
## 2 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_10pp  -0.0431 L_10pp_ROI

Catch13

catch13.modelresults <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/catch22_glasser360/svm/SVM_modelresults_catch13.Rdata")
prediction_output <- catch13.modelresults[[1]]
W <- catch13.modelresults[[2]]
SVM.weights.glasser <- as.data.frame(t(W))
SVM.weights.glasser$origlabel <- row.names(SVM.weights.glasser)
SVM.weights.glasser$label <- NA
for(weight in c(1:232)){
  newlabel <- substr(SVM.weights.glasser$origlabel[weight],1,nchar(SVM.weights.glasser$origlabel[weight])-4)
  if(substring(newlabel, 1, 1) == "R"){
    newlabelhemi <- gsub("R_", "rh_R_", newlabel)}
  else
    newlabelhemi <- gsub("L_", "lh_L_", newlabel)
SVM.weights.glasser[weight,3] <- newlabelhemi
}
SVM.weights.glasser.SA <- merge(SVM.weights.glasser, SAaxis, by="label")
cor.test(SVM.weights.glasser.SA$V1, SVM.weights.glasser.SA$final.rank, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.glasser.SA$V1 and SVM.weights.glasser.SA$final.rank
## S = 214410, p-value = 0.2533
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1083493
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_p9.46v"] <- "rh_R_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_a9.46v"] <- "rh_R_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_9.46d"] <- "rh_R_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "rh_R_OP2.3"] <- "rh_R_OP2-3"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_p9.46v"] <- "lh_L_p9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_a9.46v"] <- "lh_L_a9-46v"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_9.46d"] <- "lh_L_9-46d"
SVM.weights.glasser$label[SVM.weights.glasser$label == "lh_L_OP2.3"] <- "lh_L_OP2-3"
ggseg(.data = SVM.weights.glasser, atlas = "glasser",  mapping=aes(fill=V1), position = c("stacked")) + theme_void() + scale_fill_gradient2(high= "goldenrod1", mid = "white", low = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL)
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  side  region label         V1 origlabel 
##   <chr> <chr> <chr> <chr> <chr>  <chr>      <dbl> <chr>     
## 1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_10pp -0.499 L_10pp_ROI

Schaefer

Drug v. Placebo SVM (Schaefer)

regionheaders <- schaefer.parcel.labels$label
demoheaders <- c("subid","sesid","drug","meanFD")
colheaders <- as.matrix(c(demoheaders,regionheaders))

df <- REHO.schaefer
colnames(df) <- colheaders
df <- df %>% select(-meanFD)
df <- df[, unlist(lapply(df, function(x) !all(is.na(x))))]
model_results <- SVM_2class(df, "LOO")
saveRDS(model_results, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.schaefer.rds")
model_results <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.modelresults.schaefer.rds")
prediction_output <- model_results[[1]]
W<-model_results[[2]]
svm.model <- model_results[[3]]
num_features <- sum(!is.na(W[1,]))

Model accuracy

accuracy_fun <- function(x) sum(x$model.pred==x$drug)/dim(x)[1] # function to calculate the accuracy. 
accuracies <- sapply(prediction_output,accuracy_fun) #Apply the function to each repetition of the cross-validation.
accuracy <- mean(accuracies)
accuracy
## [1] 0.6309524

Number of correct classifications

num_obs <- length(prediction_output[[1]]$model.pred)
num_correct <- round(accuracy*num_obs)
num_correct
## [1] 53

Binomial p-value significance

b <- binom.test(num_correct,num_obs,.5) #This is a binomial test. The p-value is not used (permutation test is used instead), but this function summarizes the data nicely.
b
## 
##  Exact binomial test
## 
## data:  num_correct and num_obs
## number of successes = 53, number of trials = 84, p-value = 0.02138
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5186641 0.7337192
## sample estimates:
## probability of success 
##              0.6309524
b$pred_data <- prediction_output
b$accuracy <- accuracy

Permutation p-value significance

num_permutations = 1000
nw=1
folds="LOO"

perm_acc <- matrix(nrow = num_permutations)
perm_W <- matrix(nrow = num_permutations,ncol = num_features)
perm_list <- list()

perm_list = foreach(perm_chunk = idiv(num_permutations,chunks = nw),
                        .combine=c,
                        .export = c("featureExtraction","SVM_2class"),
                        .packages = c("dplyr","e1071","matrixTests")) %do% {  # This must be `dopar` to be parallel.
                          perm_result_list=list()
                          for (p in 1:perm_chunk){
                            # thisPerm <- perms[p,]
                            cat(sprintf("\nPermutation %d",p))
                            thisDF <- df
                            # thisDF$drug <- df$drug[thisPerm]
                            permuted <- df %>% dplyr::select(subid,drug) %>% group_by(subid) %>% mutate(perm_drug=sample(drug))
                            thisDF$drug <- permuted$perm_drug
                            perm_pred_result <- SVM_2class(thisDF,folds,
                                                           feature_selection = feature_selection,
                                                           feature_proportion = feature_proportion,
                                                           num_repetitions = 1)
                        
                            perm_result_list[[p]] = list(perm_pred_result)
                          }
                          perm_result_list
                        }

  #organizing the output
pred_data_list <- lapply(perm_list, function(x) x[[1]][[1]])
perm_Ws <- lapply(perm_list,function(x) x[[1]][[2]])
# Calculating permutation accuracies and AUCs
perm_acc_distribution <- sapply(pred_data_list, function(x) sum(x[[1]]$model.pred==x[[1]]$drug)/length(x[[1]]$drug))

perm_p <- sum(perm_acc_distribution>accuracy)/length(perm_acc_distribution)

perm_auc_distribution <- sapply(pred_data_list, function(x) AUC(DecisionValues=x[[1]]$decisionValues,labels=x[[1]]$drug))
    
    
W_test <- sapply(perm_Ws, FUN = function(x) {
    abs(x)>abs(W)
    })
W_sig <- rowMeans(W_test)
    
b$perm_p <- perm_p
b$perm_W_sig <- W_sig
b$perm_accs = perm_acc_distribution
b$perm_aucs = perm_auc_distribution


saveRDS(b, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.schaefer.rds")

saveRDS(W, file="/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.schaefer.rds")
b <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.permutations.schaefer.rds")
W <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/Inhibition/Alpraz/REHO/reho.SVM.weights.schaefer.rds")

print("How many times is the permuted accuracy greater than the actual accuracy?")
## [1] "How many times is the permuted accuracy greater than the actual accuracy?"
sum(b$perm_accs>accuracy)
## [1] 2
print(sprintf("Overall Accuracy: %1.3f; p = %.5f\n\n",accuracy,b$perm_p))
## [1] "Overall Accuracy: 0.631; p = 0.00200\n\n"

P-value distribution

hist(b$perm_accs, col="grey")
abline(v = accuracy, col="dodgerblue3", lwd=3, lty=2)

Model weights S-A axis correlation

SVM.weights.schaefer <- as.data.frame(t(W))
SVM.weights.schaefer$label <- row.names(SVM.weights.schaefer)
mappys$label <- schaefer.parcel.labels$label

SVM.weights.schaefer.SA <- merge(SVM.weights.schaefer, mappys, by="label")
cor.test(SVM.weights.schaefer.SA$V1, SVM.weights.schaefer.SA$SA_1, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$SA_1
## S = 2737344, p-value = 2.493e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4041435

Model weights inhibitory interneuron correlation

cor.test(SVM.weights.schaefer.SA$V1, SVM.weights.schaefer.SA$PVALB, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$PVALB
## S = 1504824, p-value = 0.0005336
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.228088
cor.test(SVM.weights.schaefer.SA$V1, SVM.weights.schaefer.SA$SST, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$SST
## S = 2286410, p-value = 0.009073
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1728331
cor.test(SVM.weights.schaefer.SA$V1, SVM.weights.schaefer.SA$GABRA3, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$GABRA3
## S = 2280092, p-value = 0.01048
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1695922

Model weights NMDA receptor subtype distribution

cor.test(SVM.weights.schaefer.SA$V1, SVM.weights.schaefer.SA$GRIN2B, method=c("spearman"), exact=F)
## 
##  Spearman's rank correlation rho
## 
## data:  SVM.weights.schaefer.SA$V1 and SVM.weights.schaefer.SA$GRIN2B
## S = 2187336, p-value = 0.0665
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1220123